%%
clear
L=3000;
N=512;
df=120000;
T=1/df;%常规符号持续时间
T_bar=1/df*5/4;
M=1/T_bar*1e-3;
Ts=1/N/df;
fe=-3000:50:3000;%理论
fe1=500;%仿真
te=-10/df/N:6/df/50000:10/df/N;
te1=-5/df/N:6/df/9000:5/df/N;
d=0;
for ll=1:M
    d=d+exp(1j*2*pi*ll*fe1*T_bar);
end
dt=0;
dff=1/(M/96*1e-3)/2;%相关器间隔
N0=228;
for ll=1:N0
    dt=dt+exp(-1j*2*pi*ll*te1*df);
end
fe0=-1500:10:1500;%载波频偏

fel=fe0-dff/2;%迟路
feu=fe0+dff/2;%早路
d0=abs(sinc(M*fe*T_bar));
dl=abs(sinc(M*fel*T_bar));%迟路相关值
du=abs(sinc(M*feu*T_bar));%早路相关值
Q1=dl./du-(dff-2*fe0)./(dff+2*fe0);
% Q1(ff)=du-dl;
% Q2(ff)=du+dl;
Q3=(du-dl)./(du+dl)*dff/2;
% plot(fe0,Q2)
% figure
% plot(fe0,Q1)
hold on
plot(fe,d0,'-')

xlabel('f_e (Hz)');
ylabel('f_e估计值 (Hz)');
% xlabel('f_e (Hz)');
% ylabel('f_e估计值 (Hz)');
%%
d1=sinc(N0*te*df);
figure
plot(fe,abs(d0))
hold on
% plot(fe1,abs(d)/max(abs(d)),'o')
plot(fel,abs(dl),'o')
plot(feu,abs(du),'o')
xlabel('f_e (Hz)');
ylabel('Magnitude');
grid on
% hold on
% plot(te*df*N,abs(d1))
% hold on
% plot(te1*df*N,abs(dt)/max(abs(dt)),'o')
% xlabel('f_e (Hz)');
% ylabel('Magnitude');
% grid on
%%
% 定义X和Y的范围
[X, Y] = meshgrid(fe, te);
% 定义Z为X和Y的函数，例如Z = sin(sqrt(X^2 + Y^2))
Z =(abs(sinc(M*X*T)).*abs(sinc(N0*Y*df)));
% 绘制三维曲面图
figure
surf(X, Y*df*N, Z);
% 添加标题和标签
xlabel('f_e (Hz)');
ylabel('\tau_e (T_s)');
zlabel('Magnitude');
box on
% 调整视角
% 添加颜色
colorbar; % 显示颜色条
%%
n = 100;  % 假设我们有100个迭代步骤
h = waitbar(0, 'Processing...');  % 初始化进度条

for i = 1:n
    pause(0.05);  % 模拟一些耗时操作
    waitbar(i/n, h, sprintf('Progress: %d%%', round(i/n*100)));  % 更新进度条
end

close(h);  % 关闭进度条
%% 频率RMSE验证
clear
N=512;
df=120000;
Tcoh=2e-3;%ms
T_bar=1/df*5/4;
M=Tcoh/T_bar;
df=120000;%子载波间隔（Hz）
N0=228;
CN0dB_all = 30:2.5:70;    %载噪比dbHz

for ll=1:length(CN0dB_all)
    CN0dB = CN0dB_all;
    CN0 = 10.^(CN0dB/10);
    omeg0 = 1;%相关累加前噪声方差
    C = sqrt(2*CN0/(2*N*df));%根号下相干积分前信号能量
    omeg = 1/(M*N0);%相关累加后噪声方差
    fe=-2/Tcoh/2:0.05:2/Tcoh/2;%Hz
    for ii=1:100
        NR=sqrt(omeg)*randn(1,length(fe));
        NB=sqrt(omeg)*randn(1,length(fe));
        DD=C(ll)*sinc(fe*Tcoh)+NR;
        U=abs(C(ll)*sinc(0.5+fe*Tcoh)+NR)-abs(C(ll)*sinc(0.5-fe*Tcoh)+NB);
        L=abs(C(ll)*sinc(0.5+fe*Tcoh)+NR)+abs(C(ll)*sinc(0.5-fe*Tcoh)+NB);
        D=1/(2*Tcoh)*U./L;
        fe=fe+D;
        %D2=fe-pi*(NR-NB)./(16*Tcoh*C(ll)).*(1-4*fe.^2*Tcoh.^2);    
    end
    De=D+fe;
    %De2=D2-fe;
    V(ll)=std(fe);
end
omeg1=omeg0/(N*df);%噪声功率谱密度
sigma1=sqrt(pi^2./(Tcoh^2*32*C.^2).*omeg*N0/N);
sigma1=sqrt(pi^2./(Tcoh^2*32*C.^2).*omeg);
sigma2=sqrt(pi^2*N*df*omeg1./(Tcoh^2*32*C.^2*M*N0));
%%
hold on
semilogy(CN0dB_all,V,'ro')
%%
hold on
semilogy(CN0dB_all,sigma1,'black')    
xlabel('CNR（dBHz）')
ylabel('RMSE_{f_e}（Hz）')
%% 频率二维影响
clear
L=3000;
N=512;
df=120000;
T=1/df;%常规符号持续时间
M_all=48:1:96*5;
T_bar=1/df*5/4;
Tcoh_all=M_all*T_bar;
Ts=1/N/df;
df=120000;%子载波间隔（Hz）
N0_all=52:512;
CN0dB_all = 43:2.5:43;    %载噪比dbHz

for ll=1:length(M_all)
    Tcoh=Tcoh_all(ll);
    M=M_all(ll);
    for nn=1:length(N0_all)
        N0=N0_all(nn);
        CN0dB = CN0dB_all;
        CN0 = 10.^(CN0dB/10);
        omeg0 = 1;%相关累加前噪声方差
        C = sqrt(2*CN0/(2*N*df));%根号下相干积分前信号能量
        omeg = 1/(M*N0);%相关累加后噪声方差
        fe=-1/Tcoh/2:1/Tcoh/100:1/Tcoh/2;%Hz
        for ii=1:10
            NR=sqrt(omeg)*randn(1,length(fe));
            NB=sqrt(omeg)*randn(1,length(fe));
            U=abs(C*sinc(0.5+fe*Tcoh)+NR)-abs(C*sinc(0.5-fe*Tcoh)+NB);
            L=abs(C*sinc(0.5+fe*Tcoh)+NR)+abs(C*sinc(0.5-fe*Tcoh)+NB);
            D=1/(2*Tcoh)*U./L;
            fe=fe+D;
            %D2=fe-pi*(NR-NB)./(16*Tcoh*C(ll)).*(1-4*fe.^2*Tcoh.^2);    
        end
        De=D+fe;
        %De2=D2-fe;
        V(ll,nn)=std(fe);
    end
end
[X, Y] = meshgrid(M_all, N0_all);
% 定义Z为X和Y的函数，例如Z = sin(sqrt(X^2 + Y^2))
Z=V';
% 绘制三维曲面图
figure
surf(X/480, Y/512, Z/max(max(Z)));
% 添加标题和标签
xlabel('\alpha_\tau');
ylabel('\alpha_f');
zlabel('RMSE_{f_e}(Hz)');
% set(gca,'ZLim',[0 1400]);%X轴的数据显示范围
% set(gca,'ZTick',[0:100:1400]);%设置要显示坐标刻度
hold on
Z=Z./Z/max(max(Z))*100;
surf(X/480, Y/512, Z);
box on

%%
omeg1=omeg0/(N*df);%噪声功率谱密度
sigma1=sqrt(pi^2./(Tcoh^2*32*C.^2).*omeg);
sigma2=sqrt(pi^2*N*df*omeg1./(Tcoh^2*32*C.^2*M*N0));
hold on
semilogy(CN0dB_all,V,'ro')
hold on
semilogy(CN0dB_all,sigma1,'black')    
xlabel('CNR（dBHz）')
ylabel('RMSE_{f_e}（Hz）')

%%
figure
plot(fe,-D)
%% ifft
clear
a=zeros(1,100);
a(1:3)=1;
b=ifft(a);
%%
%% 时延环路
clear
L=3000;
N=512;
df=120000;
T=1/df;%常规符号持续时间
T_bar=1/df*5/4;
M=1/T_bar*1e-3;
Ts=1/N/df;
fe=-3000:50:3000;%理论
fe1=500;%仿真
te=-10/df/N:6/df/50000:10/df/N;
te1=-5/df/N:6/df/9000:5/df/N;
d=0;
N0=228;
dt=0;
dff=1/(N0*df)*1.5;%相关器间隔
N0=228;
for ll=1:N0
    dt=dt+exp(-1j*2*pi*ll*te1*df);
end
te0=-2/df/N:6/df/300000:2/df/N;%载波频偏

del=te0-dff/2;%迟路
deu=te0+dff/2;%早路
d0=abs(sinc(N0*te*df));
dl=abs(sinc(N0*del*df));%迟路相关值
du=abs(sinc(N0*deu*df));%早路相关值
%Q1=dl./du-(dff-2*fe0)./(dff+2*fe0);
% Q1(ff)=du-dl;
% Q2(ff)=du+dl;
Q3=(du-dl)./(du+dl)*dff/2/Ts;
% plot(fe0,Q2)
% figure
% plot(fe0,Q1)
hold on
plot(te0/Ts,-Q3,'-')
% figure
% plot(te0/Ts,du,'-')
xlabel('\tau_e (Hz)');
ylabel('\tau_e估计值 (Hz)');
%% 时延RMSE验证
clear
N=512;
df=120000;
Tcoh=5e-3;%ms
T_bar=1/df*5/4;
M=Tcoh/T_bar;
N0=228;
CN0dB_all = 40:3:70;    %载噪比dbHz
dtt=1/(N0*df);%相关器间隔
for ll=1:length(CN0dB_all)
    CN0dB = CN0dB_all;
    CN0 = 10.^(CN0dB/10);
    omeg0 = 1;%相关累加前噪声方差
    C = sqrt(2*CN0/(2*N*df));%根号下相干积分前信号能量
    omeg = 1/(M*N0);%相关累加后噪声方差
    te0=-1/df/N:6/df/3000000:1/df/N;%时延偏差
    for ii=1:10
        del=te0-dtt/2;%迟路
        deu=te0+dtt/2;%早路
        NR=sqrt(omeg)*randn(1,length(te0));
        NB=sqrt(omeg)*randn(1,length(te0));
        U=abs(sqrt(N/N0)*C(ll)*sinc(N0*del*df)+NR)-abs(sqrt(N/N0)*C(ll)*sinc(N0*deu*df)+NB);
        L=abs(sqrt(N/N0)*C(ll)*sinc(N0*del*df)+NR)+abs(sqrt(N/N0)*C(ll)*sinc(N0*deu*df)+NB);
        % U=abs(C(ll)*sinc(N0*del*df)+NR)-abs(C(ll)*sinc(N0*deu*df)+NB);
        % L=abs(C(ll)*sinc(N0*del*df)+NR)+abs(C(ll)*sinc(N0*deu*df)+NB);
        D=dtt/(2)*U./L;
        te0=te0-D;
    end
    De=D+te0;
    %De2=D2-fe;
    V(ll)=std(te0);
end
omeg1=omeg0/(N*df);%噪声功率谱密度
sigma1=sqrt(pi^2./((N0*df)^2*32*C.^2).*omeg*N0/N);
% sigma1=sqrt(pi^2./((N0*df)^2*32*C.^2).*omeg);
%sigma2=sqrt(pi^2*N*df*omeg1./(Bcoh^2*32*C.^2*M*N0));
hold on
semilogy(CN0dB_all,V*3*1e8,'r-o')
%semilogy(CN0dB_all,sigma1*3*1e8,'black')    
xlabel('CNR（dBHz）')
ylabel('RMSE_{t_e}（m）')

%% 时延二维影响
clear
L=3000;
N=512;
df=120000;
T=1/df;%常规符号持续时间
M_all=48:1:96*5;
T_bar=1/df*5/4;
Tcoh_all=M_all*T_bar;
Ts=1/N/df;
df=120000;%子载波间隔（Hz）
N0_all=52:512;
CN0dB_all = 43:2.5:43;    %载噪比dbHz
for ll=1:length(M_all)
    Tcoh=Tcoh_all(ll);
    M=M_all(ll);
    for nn=1:length(N0_all)
        N0=N0_all(nn);
        dtt=1/(N0*df);%相关器间隔
        CN0dB = CN0dB_all;
        CN0 = 10.^(CN0dB/10);
        omeg0 = 1;%相关累加前噪声方差
        C = sqrt(2*CN0/(2*N*df));%根号下相干积分前信号能量
        omeg = 1/(M*N0);%相关累加后噪声方差
        te0=-1/df/N:2/df/N/100:1/df/N;%时延偏差
        for ii=1:10
            del=te0-dtt/2;%迟路
            deu=te0+dtt/2;%早路
            NR=sqrt(omeg)*randn(1,length(te0));
            NB=sqrt(omeg)*randn(1,length(te0));
            U=abs(sqrt(N/N0)*C*sinc(N0*del*df)+NR)-abs(sqrt(N/N0)*C*sinc(N0*deu*df)+NB);
            L=abs(sqrt(N/N0)*C*sinc(N0*del*df)+NR)+abs(sqrt(N/N0)*C*sinc(N0*deu*df)+NB);
            U=abs(C*sinc(N0*del*df)+NR)-abs(C*sinc(N0*deu*df)+NB);
            L=abs(C*sinc(N0*del*df)+NR)+abs(C*sinc(N0*deu*df)+NB);
            D=dtt/(2)*U./L;
            te0=te0-D;
        end
        %De2=D2-fe;
        V(ll,nn)=std(te0);
    end
end
[X, Y] = meshgrid(M_all, N0_all);
% 定义Z为X和Y的函数，例如Z = sin(sqrt(X^2 + Y^2))
Z=V';
% 绘制三维曲面图

hold on
Z=Z*3e8;
surf(X/480, Y/512, Z/max(max(Z))+1);
% 添加标题和标签
xlabel('\alpha_\tau');
ylabel('\alpha_f');
zlabel('RMSE_{\tau_e}(m)');
%colormap(cool);   % 更改色图为热图色系
set(gca,'XLim',[0.1 1]);%X轴的数据显示范围
%set(gca,'ZTick',[0:100:1400]);%设置要显示坐标刻度
hold on
Z=Z./Z/max(max(Z))*3+1;
surf(X/480, Y/512, Z);
box on
%%
clear
L=3000;
N=512;
df=120000;
T=1/df;%常规符号持续时间
M_all=48:1:96*5;
T_bar=1/df*5/4;
Tcoh_all=M_all*T_bar;
Ts=1/N/df;
df=120000;%子载波间隔（Hz）
N0_all=52:512;
CN0dB_all = 43:2.5:43;    %载噪比dbHz
CN0 = 10.^(CN0dB_all/10);
lam_f=50;%Hz
lam_t=1;%m
lam_t=lam_t/3/1e8;%s
M0=1;
N0=1;
for ll=1:800
    M_all(ll)=M0;
    N_all(ll)=N0;
    sigmat=sqrt(pi^2*N*df/(32*CN0*M0*N0*(N0*df)^2));
    sigmaf=sqrt(pi^2*N*df/(32*CN0*M0*N0*(M0*T_bar)^2));
    if sigmat>lam_t
        N0=N0+1;
    else
        N0=N0-1;
        M0=M0-1;
    end
    if sigmaf>lam_f
        M0=M0+1;
    else
        M0=M0-1;
        N0=N0-1;
    end
end
plot(M_all);
hold on
plot(N_all);
%%
M0=sqrt(pi/(T_bar^1.5*sqrt(32*CN0/(N*df*df))*sqrt(lam_f^3/lam_t)));
N0=sqrt(pi/(df^1.5*sqrt(32*CN0/(T_bar*N*df))*sqrt(lam_t^3/lam_f)));
M0=ones(1,800)*M0;
N0=ones(1,800)*N0;
hold on
plot(M0);
hold on
plot(N0);
%%
%% 频率RMSE验证
clear
N=512;
df=120000;
Tcoh=1e-3;%ms
T_bar=1/df*5/4;
M=Tcoh/T_bar;
df=120000;%子载波间隔（Hz）
N0=228;
CN0dB_all = 60:2.5:60;    %载噪比dbHz
mentor=100;
for ll=1:mentor
    CN0dB = CN0dB_all;
    CN0 = 10.^(CN0dB/10);
    omeg0 = 1;%相关累加前噪声方差
    C = sqrt(2*CN0/(2*N*df));%根号下相干积分前信号能量
    omeg = 1/(M*N0);%相关累加后噪声方差
    fe=3/Tcoh;%Hz
    for ii=1:1000
        fe_all(ll,ii)=fe;
        NR=sqrt(omeg)*randn(1,length(fe));
        NB=sqrt(omeg)*randn(1,length(fe));
        DD=C*sinc(fe*Tcoh)+NR;
        dff=0.6/Tcoh;
        U=abs(C*sinc(0.5*Tcoh*dff+fe*Tcoh)+NR)-abs(C*sinc(0.5*Tcoh*dff-fe*Tcoh)+NB);
        L=abs(C*sinc(0.5*Tcoh*dff+fe*Tcoh)+NR)+abs(C*sinc(0.5*Tcoh*dff-fe*Tcoh)+NB);
        D=1/(2*Tcoh)*U./L;
        fe=fe+D;
        %D2=fe-pi*(NR-NB)./(16*Tcoh*C(ll)).*(1-4*fe.^2*Tcoh.^2);    
    end
    % De=D+fe;
    % %De2=D2-fe;
    % V(ll)=std(fe);
end
fe_all=mean(fe_all);
hold on
plot(fe_all);
% omeg1=omeg0/(N*df);%噪声功率谱密度
% sigma1=sqrt(pi^2./(Tcoh^2*32*C.^2).*omeg*N0/N);
% sigma1=sqrt(pi^2./(Tcoh^2*32*C.^2).*omeg);
% sigma2=sqrt(pi^2*N*df*omeg1./(Tcoh^2*32*C.^2*M*N0));
% figure
% semilogy(CN0dB_all,V,'ro')
% hold on
% semilogy(CN0dB_all,sigma1,'black')    
% xlabel('CNR（dBHz）')
% ylabel('RMSE_{f_e}（Hz）')
%%
x=-2:0.1:2;
y=x;
hold on
plot(x,y);